Classical density-functional theory for water 
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We introduce a new computationally efficient and accurate classical density-functional theory for 
water and apply it to hydration of hard spheres and inert gas atoms. We find good agreement 
with molecular dynamics simulations for the hydration of hard spheres and promising agreement for 
the solvation of inert gas atoms in water. Finally, we explore the importance of the orientational 
ambiguity in state-of-the-art continuum theories of water, which are based on the molecular density 
only. 



INTRODUCTION 

Water is the most abundant liquid on earth and is es- 
sential for life. Water is a called a universal solvent for its 
ability to dissolve both acids and bases, as well as polar 
and non-polar molecules. 

Despite its importance to a wide range of problems 
in chemistry and biology, it is still an open challenge to 
build a first principles theory of the microscopic structure 
of water and its interaction with solutes [if. The underly- 
ing reason for this is the complex interplay between the 
kinetic energy of molecules and the the orientation de- 
pendent potential energy of inter-molecule interactions. 

For the study of liquids and their interactions with ex- 
ternal systems, two general classes of theoretical meth- 
ods have been developed. The first treats the liquid as 
a collection of molecules treated either ab initio within 
density-functional theory or with a model interatomic po- 
tential, such as the simple point charge (SPC) model [H 
or the 4-site transferable intermolecular potential TIP4P 
and then uses molecular dynamics (MD) or Monte 
Carlo numerical methods to perform statistical averag- 
ing over the thermal phase space These methods 
are intuitively simple to understand and also relatively 
straightforward to implement numerically. However, be- 
cause they involve statistical averaging of many molecules 
over an exponentially growing configuration phase space, 
these methods arc numerically very demanding. The sec- 
ond class of methods treats the liquid as a continuous 
media [7- 11 [. Without the need for thermal averaging 
or to represent molecules individually, these latter meth- 
ods are much more efficient computationally. However, 
such models are generally built in an empirical way and, 
generally, there is no systematic way to improve them. 
Nonetheless, significant progress has been made in un- 
derstanding the interaction of water with external envi- 
ronments using this approach, such as the work of Pratt 
and Chandler 12[ on the theory of the hydrophobic ef- 
fect and the Lum, Chandler, Weeks (LCW) theory of 
hvdrophobicitv [1 3j . 

In this work, we explore the nature of the interaction 
of water with external environments using a somewhat 
different approach than that of Chandler and coworkers 



and work in a classical density-functional theory frame- 
work. There are two main advantages of working in this 
framework. First, the classical density- functional the- 
ory of liquids is a continuous theory of the liquid state 
which is exact in principle. Moreover, this framework 
gives the free energy and the density profile of the liq- 
uid in any external potential V(r) in terms of a single 
density- functional [lj , so that study of the hydrophobic 
effect (liquid in contact with either an impenetrable wall 
or an impenetrable hard sphere) and of the interaction of 
the liquid with any solute can be carried out in a single, 
unified framework. A number of approximate density- 
functionals have been developed for water and applied to 
the hydrophobic effect [Tol hj. Much more demanding 
theories which go beyond the average molecular density 
to consider explicitly the distribution of molecular orien- 
tations in water have also been developed HI 1(| . All of 
these density functional theories, however, prove to be 
either quite computationally demanding or provide an 
over-simplified description. 

We begin this work by introducing a new, compu- 
tationally efficient density-functional theory for water 
which accurately reproduces the hydrophobic effect near 
hard boundaries. We then present the first application 
of a classical density-functional theory to realistic ab ini- 
tio potential energy surfaces of solutes, applying our the- 
ory to the solvation of the inert gas sequence. This lat- 
ter study allows us to address directly the question of 
whether explicit orientation dependence, with all of the 
concomitant computational demands, is necessary to pro- 
vide an accurate description of the solvation of even the 
simplest solutes. 



CLASSICAL DENSITY-FUNCTIONAL THEORY 

Classical density-functional theory for liquids follows 
from a variational principle which was first established by 
Hohenberg and Kohn[17| for the inhomogeneous electron 
gas at zero temperature and later generalized to finite 
temperature by Merminflij. 

The theorem states that the equilibrium grand poten- 
tial fio of a liquid at chemical potential fi in any exter- 
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nal potential V(r) can be found by minimizing, over all 
possible thermally averaged molecular densities p(r), the 
sum of a universal functional and the interaction of that 
density with the potential, 



fio = min ^ F[p(r)] + / (V(f) — p) p(r) dr 

p(r' 



(1) 



with the density profile p(r) which minimizes the free- 
energy giving the exact average density at thermal equi- 
librium. F[p(r)) is said to be a universal functional of 
density p(r) in the sense that it does not depend what- 
soever on the external potential V(r) but only upon the 
nature of the solvent itself. This yields the significant 
advantage that once a good approximation is found for 
the universal functional _F[p(r)], it can be used to study 
the interaction with any possible external potential V(r), 
such as that representing one solute or another. 

Even though F[p(f)] is exact in principle, its analytical 
form is not known and approximations must be devel- 
oped. Generally, one begins by separating the functional 
as 



F[p(r)] = F id [p(f)}+F ex [p(r)}, 



(2) 



where Fi d [p(r)] is the free energy of ideal gas whose an- 
alytic form is known exactly, 



F ld [p(f)} = kT J p(f)\n(p(f))df, 



(3) 



and where F ex [p(r)\ is defined as the excess free energy 
beyond that of the ideal gas and contains all the other 
terms responsible for correlations among the molecules. 

Our density-functional theory is inspired by a class of 
weighted-density functionalsQ- Our method, however, is 
computationally much simpler in that it does not require 
computationally demanding self-consistent calculations 
of the weighted density. On the other hand, our form 
does allow us to incorporate much of the same physics 
as [3] and thus find a form which is competitive compu- 
tationally but more accurate than other, more simplified 
functional [10|. 



CONSTRUCTION OF CLASSICAL 
DENSITY-FUNCTIONAL THEORY FOR WATER 

The form we choose to consider is 

F[p(r)} = F ld [p(r)} + J fe*Q>(f))dF (4) 

+ ff MxP- 

where Fi d [p(r)] is given by Eq. and f ex (p(r)) is a local 
function of p(r), a Gaussian smoothed density 



Vapor density 


0.023 kg/m s 


Water density 


997.1 kg/m 3 


Bulk modulus 


2.187 GPa 


dB/dP 


5.83 



6.630 x 10 J 



-1.277 x 10 ; 



9.200 x 10 5 



-2.602 x 10 J 



e 


8.906 x 10" 4 


f 


-1.415 x 10~ 2 


g 


1.077 x 10~ 10 



p(f) = (2Yn)-i / p{f)e df , (5) 



TABLE I: Properties of uniform bulk water at Standard Am- 
bient Temperature and Pressure (SATP: T=25 °C, P=100.00 
kPa) 



and x(f) is an interaction kernel which integrates to zero. 
The first term in Eq. (j4]) , which is known exactly, ensures 
that the short length-scale properties (to which F ex does 
not couple) are treated properly. The second term in 
Eq. ((4]) ensures that a large-scale phase separation be- 
tween vapor and liquid is treated correctly and, often, 
can capture exactly the experimental two-particle corre- 
lation functions in the bulk liquid [7] for an appropriate 
convolution kernel in Eq. ([5]). However, we find that, 
for the experimentally observed two-particle correlation 
function in liquid water, the equation for the kernel has 
no real solutions. The solutions have small imaginary 
parts and thus cannot be used directly. The real parts of 
the solution, however, closely resemble a Gaussian. Thus, 
we employ a Gaussian in Eq. ([5]) and add the last term 
in Eq. (j4|) to ensure that the experimental two-particle 
correlation function in the bulk liquid is reproduced. 

To construct the functions and parameters needed in 
Eq. Q, we begin by noting that, for the uniform liquid, 
the last term gives zero (x has zero integral) and p(r) — 
p{r), so that f e x(p) can be constructed to reproduce the 
properties of the uniform liquid exactly. In practice, we 
parameterized f e x(p) as a sixth-order polynomial 

f ex [p] = ap 6 + bp 5 + cp 4 + dp 3 + ep 2 + fp + g. (6) 

with parameters chosen to reproduce, for bulk water at 
standard ambient temperature and pressure (SATP), the 
density and bulk modulus of both the liquid and vapor 
(four parameters), the derivative of the bulk modulus of 
the liquid with respect to pressure dB/dP, and phase 
coexistence of the liquid and vapor (equality of grand 
potential) . Table Q] gives both the fit data and the re- 
sulting parameters, and Figure Q] presents the final func- 
tion F(p). Note that in this work we employ mostly 
atomic units (a.u.) so that Planck's constant, the elec- 
tron mass and the fundamental charge all have numerical 
value unity (h = m e = e = 1), implying that energies are 
measured in units of 1 hartree=27.21 eV and distances 
in units of 1 bohr=0.5291 A. 

To construct x( 7? )i we work in the Fourier represen- 
tation. As noted above, %(r) integrates to zero so that 
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FIG. 1: Free energy of water as a function of its density at 
SATP 



x{k = 0) = 0. For k ^ 0, we fit to the experimen- 
tal pair distribution function g(r) of bulk uniform wa- 
ter at SATP[19(. Finally, the Gaussian width parameter 
A = 0.4415 a.u. has been fit to the surface tension of bulk 
water 70 mN/m or 4.5 x 10~ 5 a.u. We emphasize that 
in constructing this functional we used only the macro- 
scopic properties of water. 



APPLICATION TO HYDRATION OF HARD 
SPHERES 



To test the accuracy of this approximation we use it to 
calculate the free energy of a spherical cavity in water, a 

" to 



standard test case used in the literature E, E3, S EH 



explore density variations in water over all possible length 
scales. The interaction potential with a hard sphere is 
over-idealized because it depends only on the distance 
of the water molecule to the hard sphere and does not 
depend on the orientation. For instance, for a real water 
molecule it is unclear even what point to take to represent 
the location of the molecule. The position of the oxygen 
nucleus is most often taken. 

Figure [5] compares our results for the surface ten- 
sion (free energy change per unit area) of a spherical 
cavity with results of molecular simulations for SPC/E 
water [2^. The figure verifies that the surface tension ap- 
proaches in the macroscopic value in the limit of large 
radii and show the, for smaller radii, the surface ten- 
sion has a strong dependence on radius. Our results are 
in agreement with Lum- Chandler- Weeks [l^ theory. The 
good agreement with explicit molecular dynamics simu- 
lations suggests that this model gives good quantitative 
description of the hydrophobic effect, a central problem 
of theoretical chemistry. 
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FIG. 2: Surface tension of a hard sphere in water at 
SATP. Solid line denotes the results of our classical density- 
functional theory calculations. Molecular simulation results 
for a cavity in SPC/E water are indicated by diamonds [221 ], 
Dashed line corresponds to macroscopic surface tension of wa- 
ter. 



APPLICATION TO HYDRATION OF INERT GAS 
ATOMS 

The ultimate motivation for inhomogeneous continu- 
ous theories of bulk water is to understand solvation of 
real solutes, not artificial hard boundaries. Within the 
density- functional framework of Eq. ([1]) , one can simply 
and rigorously incorporate the effects of any external po- 
tential V(r) acting on the liquid, making treatment of 
real solutes simple in principle, provided a static poten- 
tial V(f) accurately describes the interaction of water 
with the solute. For our next test, we treat a simple but 
challenging problem for which there is experimental data, 
the solvation of inert gas atoms. 

The experimental solubilities of the inert gasses in wa- 
ter give the corresponding solvation free energies directly. 
In the dilute limit, the free energy of solvation Af2 relates 
to the mole fraction solubility through the relation 



Afl = fc s Tlog 



Xin v 



k B T log 



latm 



ks TX\n u 



(7) 



where X\ is the mole fraction solubility in water and n w 
is the number density of liquid water. 

The experimental data, which Table HI summarizes, 
show an interesting trend. Although the solvation en- 
ergies of hard spheres increase with size, the solvation 
energies of inert gas atoms tend to decrease. (Note that 
for helium, the worst possible case, we estimate quan- 
tum zero-point effects to be quite small, on the order of 
0.005 eV, so that this must be purely an effect of the 
interactions.) Thus, a simple cavity model for solvation 
of inert gas atoms is not a good approximation and one 
must use a more realistic interaction potential V(f). 
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Exp. (eV) 


V kT (eV) 


Vmin (eV) 


Helium 


0.12 


0.12 


0.09 


Neon 


0.12 


0.11 


0.07 


Argon 


0.09 


0.22 


0.17 


Krypton 


0.07 


0.25 


0.18 



R (Angstrom) 



TABLE II: Solvation energies of inert gas atoms 



bigger the discrepancy also increases suggesting that ori- 
entational ordering of water molecules around the atoms 
is playing an important role and a density-only descrip- 
tion is not enough. 



FIG. 3: Interaction of Argon and a water molecule. Each 
diamond correspond to a specific orientation. Also presented 
are the Boltzmann average potential (solid line) and the min- 
imum potential (dot-dashed line). 

We determined the potential energy of interaction be- 
tween a water molecule and an atom of each of the in- 
ert gasses in Table QT] directly through ab initio density- 
functional theory calculations within the generalized gra- 
dient approximation (GGA) [23|. For these calculations 
we take r*to be the displacement between the nucleus of 
the inert gas and the nucleus of the oxygen atom and con- 
sidered one hundred different orientations for the water 
molecule, with these orientations optimized to sample the 
angular dependence evenly according to the procedure of 
Womersley and Sloan [24|. Figure [3] summarizes the re- 
sults for argon, which was typical of all of the inert gasses 
considered in this work. 

The scatter in the interaction energy at each distance r 
shows a strong dependence on the orientation of the wa- 
ter molecule, particularly at closer distances. The ques- 
tion then immediately arises of how to couple to any con- 
tinuum theory based solely on the molecular density N(r) 
with no information about the orientations of the solute. 
As the coupling in such density-only functional theories 
takes the form J V(r)n(r)dr, a choice must be make to 
to define a unique value for the potential V(r) for each 
point f. One approach would be to assume that each 
water molecule independently assumes its most favored 
orientation for a given distance, thereby defining V(r) as 
the minimum energy envelop (V m i n ) of the data in Fig- 
ure [H Another approach would be to take the thermal 
average interaction under the assumption that the water 
molecules choose their orientations independently of each 
other, thereby defining V(r) as the Boltzmann weighted 
average (VkT) 01 the data at each distance. 

Table [IT] summarizes the results of the minimum free 
energy of our functional for both of these approaches. 
For helium and neon experimental and theoretical results 
agree quite well. Also the difference between using V m in 
and VkT is less than O.OAeV. As the size of atoms gets 



CONCLUSIONS AND FUTURE DIRECTIONS 

This work presents a simple density-functional theory 
for water which gives reasonable overall agreement with 
molecular dynamics simulation data for the solvation of 
hard spheres in water, with quite good agreement for 
smaller spheres (radius less than 2.5A). However, as we 
move to look at simple solutes, such as inert gasses, we 
find mixed results. While we are able to reproduce the 
counter-intuitive trend of decreasing energy with increas- 
ing size in the inert gas sequence when going from He to 
Ne, we fail to do so for the larger atoms. Thus, although 
solvation of hard spheres, whose interaction with water 
molecules is generally taken to involve only the location 
of the oxygen nucleus and thus contains no orientation 
dependence, may be well described by density-only the- 
ories, such theories do not necessarily describe well the 
solvation of even simple solutes such as inert gases. It 
thus appears that treatment of orientation in some form 
is needed to attain results approaching chemical accu- 
racy for realistic solutes. This work was funded by NSF 
GRANT #CHE-0113670. 
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